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Abstract 

After a brief introduction to the statistical description of data, these lecture notes 
focus on quantum field theories as they emerge from lattice models in the critical 
limit. For the simulation of these lattice models, Markov chain Monte-Carlo methods 
are widely used. We discuss the heat bath and, more modern, cluster algorithms. 
The Ising model is used as a concrete illustration of important concepts such as 
correspondence between a theory of branes and quantum field theory or the duality 
map between strong and weak couplings. The notes then discuss the inclusion of 
gauge symmetries in lattice models and, in particular, the continuum limit in which 
quantum Yang-Mills theories arise. 



Notes based on a lecture presented at the XIX Physics Graduate Days at the 
University of Heidelberg, 8th - 12th October 2007. 
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1 Statistical data analysis 
1.1 The central limit theorem 

Assume that we would like to determine a physical observable x such as a hadron mass 
or a decay constant by a numerical calculation involving statistical methods or by a direct 
experimental measurement. A perfect device would just produce x with a single measure- 
ment. In practice, such a device does not exist. A realistic device produce a value x in 
the interval [x, x -|- dx\ with probability P{x) dx, where the probability distribution P{x) 
characterises the apparatus. We will not assume that our experimental device is hampered 
by systematic errors, but we will assume that the device produces the exact value x by an 



2 



average over many measurements, i.e., 

jd.xp(x) = ,, (1) 

but, depending on P{x), a single measurement can be far from the true value. 

As an example, we consider an observable x = 3 and a crude experiment which can produce 
any value for between and 6 with equal probability: 

pf^) = / V6 fora;e[0,6] 
[0 otherwise. 

Obviously, a single measurement for x is not sufficient to reveal the true observable. The 
only thing we can do is to repeat the measurement n times and to consider the average: 



1 

y = - 

n 



Xi + . . . + Xr. 



where xi to Xn are the values obtained from each of the measurements. For the moment, 
we will assume that the measurements are independent, i.e., that the probability for finding 
a set {xi . . . Xn} of data is given by: 

P(Xi, ...,Xn) = -P(Xi) . . . P{Xn) ■ 

The crucial question is to which accuracy have we estimated the true observable x? 
The answer can be inferred from the probability distribution Q{y) for the value y: 

Qn{y) = Jf[dxiS(^y-^[xi + ... + x„]^ P{xi) ... P{xn) . (3) 

Given the proper normalisation of the single event distributions, i.e.. 



dxi P{xi) = 1 , 

using ([T]), we can easily show that the average of y coincides with the observable: 

dy y Qn{y) 

^ - - [xi + . . . + x„] ) Pixi) . . . PiXn) 



j dy j Wdxi^[xi + . . . + Xn] 6 - ^[2 

/W dXi — [Xi + . . . + Xn\ PiXi) ... P(Xn) = —Tl I dx X Pix) = X 

^■^ n n 
1=1 
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A natural measure for the error a of our estimate is provided by the second moment: 

i) = j dy [y - Qn{y) , (the variance) . (4) 



If the distribution Qn{y) peaks around the true value for our observable x and a{n) is tiny, 
it would mean that a single estimator y has high probability to fall close to x with high 
probability implying that it yields a good approximation to x. 



Let us study the moments of the distribution Q{y): 



dy Qn{y) y"^ • 



(5) 



In order to draw further conclusions, we need to restrict the classes of single event proba- 
bility distributions P{x): we will assume that its Fourier transform 



P[f3) = j dxP{x) exp{-2/3x} 



(6) 



is an analytic function of /5 at /5 = 0. As a consequence, the moments of P{x) exist and 
are given by: 



Pn 



dx Fix) X™ = 



dp' 



1/3=0 



(7) 



We will further assume that P{P) vanishes for \(3\ ±oo. This seems to be quite a weak 
constraint. I point out, however, that systems with rare but large fluctuations generically 
fail to possess higher moments. One example is stock market indices [1]. 

Our aim is to express the moments of Qn{y) in terms of the moments of P{x). For this 
purpose, we rewrite the (5-function in as 



and find 



6 [y [xi + 



Qn{y) 



da 

— exp{tay) 

ZTT 



da r- n TT f • 

— exp z ay\ \\exp<-i-Xi( , 
2tt -'--'-In 

4 = 1 



J dx P{x) exp|— z — 



da 
2^ 



exp(i a y) 



P{- 



The moments of Qn{y) are then obtained from 

d"" 



/ 2^ da- 



exp(i a y) 



n 



(9) 
(10) 
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After a series of partial integrations with respect to a (note that boundary terms vanish 
by virtue of our above assumptions), the latter equation is given by 



da 

dy /^exp(.ay)(.r^ 



da 1 
27r 71^ 



dy exp(i a y) (i) 



pn 
1^ 



■m 



1/3=0 



.(12) 



Of particular interest are the so-called cumulants Cfc[Q„] of the distribution Qniv)- These 
are defined via the generating function 



To(x) 



Ck[Qn] ■-- 



m=0 



d^ 

dx^ 



In Tq(x) 1^.= 



Note that in particular we find for the 'error' a in (jl]) 

= 92 - gi = C2[Qn] ■ 

Using Taylor's theorem and the explicit expression (1121) . we find 



Tq{x) 



PN — 

n 



Ck[Qr 



n 



k-l 



InP(O) 



(fc) 



(13) 



(14) 



(15) 



where (A;) denotes the kth derivative. Introducing the cumulants Cfc[P] of the single event 
distribution as well, i.e., 



Ck[P] 

we arrive at a very important result: 

CklQr 



InP(O) 



(fc) 



(16) 



(17) 



Note that the cumulants Ck[P] are finite numbers which characterise the single event proba- 
bility distribution. Equation (fT7|) then implies that for increasing number of measurements 
n, the higher [k > 1) cumulants of Qniv) vanish. In particular, we find that 

o-{n) = a/ C2 [Qn] = '^'^^i-^^ oc l/\/n 



n 



For the above example, we find 



Pi = - dx X 
and therefore 



P2 



dx X 



12 



C2[P] = 12-32 = 3, (19) 



cr 



(n) = ^/3/n . 
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It is well known that if Ck[G\ = for k > 2, the probability distribution G is a Gaussian. 
We therefore expect that if n is chosen sufficiently large so that we can neglect [Qn] with 
k > 2, we should be able to approximate Qniy) by a Gaussian. To support this claim 
(without mathematical rigour), we start from (|TOl) : 



da 



27r 



Qn{.y) 

and expand the logarithm with respect to a 

Qn{y) = 



exp(z a y) exp < n In 



P 



da 
2^ 



da 
2^ 



fc=o • n j 

{oo ^ 
fc=i 



where we have used -P(O) = J dxP{x) = 1 and the definition of the cumulants of P in 
( |T6l) . Using Ci[P] = = X = I/, we find: 



Qn(l/) 



2^ 



exp[ia{y-y)] exp 



+ C(aVn2 



(20) 



Note that the dominant contributions from the a integration arises from the regime where 
a < ^Jn. In this regime, the correction term is of order 



0(aV^' 



0{\ 



and will be neglected for sufficiently large n. The remaining integral can be easily per- 
formed: 



Qn{y) 



'2tio 



exp 



2a2 



n 



(21) 



which is the celebrated Gaussian distribution. This finding is called the central limit 
theorem: if the moments of probability distribution P{x) exist, the probability distribution 
for the average y can be approximated by a Gaussian for sufficiently large n given that the 
standard deviation a is properly scaled with n. 

Let us discuss this result. Figured] shows the distributions Qn{y) for n = 1, 2, 3 and n = 10. 
Already for n = 10 the distribution is well approximated by the Gaussian. 

The existence of at least the moment C2[P] of the single event distribution is crucial for 
the error reduction by repeated measurements. Let us consider a Lorentz distribution for 
the moment: 

Pl(x) = — . (22) 



b 1 



{x/by 



6 



0.8 




Figure 1: Illustration of the central limit theorem: probability distributions of the average 
y after n measurements. 

With the naked eye this distribution resembles the Gaussian. The crucial difference is, 
however, that the second moment does not exist: 



dx Pl ix) 



oo . 



The Fourier transform of Pl{x) does, however, exist 

PlW) = exp|-6|/3|} . (23) 

If we now repeat the measurements n times, the distribution of the average y is, according 
to (dU]), by 



Qn{y) 



da 
2^ 



exp(i a y) 



exp ( —h 



a 



n 



Obviously, the probability distribution does not change at all even if we repeat the mea- 
surements many times. This actually implies that it impossible to experimentally gain a 
reliable value for the observable x. 



1.2 Error analysis 

Let us return to the example in ([2]), and let us assume a group of experimentalists has 
performed n = 12 measurements with the result: 
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2.813 2.021 0.331 0.865 
5.394 5.937 5.027 1.556 
0.325 2.750 1.283 3.890 



The average over these values and an estimate {x^) of C2[P] are given by 

^ n 1 " 

y = -y^Xk ^ 2.683, (x^) = - V x? ^ 3.581. (24) 
n ^-^ n ^-^ 

We point out that 3.581 is a poor estimate of the true value (fT9l) of C2[P] = 3, but it gives 
the order of magnitude. With this estimate for C2[P], we find for the error ( ITSl) 

/3.581 

a(n=12) ^ ~ 0-546. 

Hence, the final 'experimental' result for the observable would be 

X ^ 2.683 ±0.546 = 2.7(5). (25) 
Note that the true result x = 3 lies well with the reach of the error bars. 

The above experiment was repeated by several research labs. Depending on the budget 
and the focus of research, different labs produce different numbers n of measurements: 



CERN GSI DESY BNL 

n 120 50 78 150 

y 3.112 ±0.163 2.764 ± 0.255 3.110 ±0.207 3.083 ± 0.143 



The smallest error was produced by the largest experiment (BNL). We could just quote 
their result, but it would be a pity to disregard a total of 248 measurements which were 
carried out by the other groups. How can we obtain a 'world average'for the observable x 
and how can we quantify its (statistical) error? 

To answer these questions, we assume that the number n of each measurement was large 
enough to approximate the distribution of an individual result yk, k = 1 . . . N (where 
= 4 for the above example) by a Gaussian fl^Tl) : 

Qiyi) - ^expj-^^^l . (26) 
For the world average y we make the ansatz 

TV TV 

y = ^aiVi , X^^i = 1 • (27) 
1=1 1=1 
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The weights ai must be chosen in an optimal way. This choice will depend on the errors ai 
of the individual experiments. In particular, the experiment with the smallest error should 
contribute to the world average with the largest weight. Assuming that the experiments at 
the different labs were carried out independently, the probability distribution of the world 
average is now given by 

„ Af N 

wiy) = lldy,6{y ~ J2^^y^) Qiyi)---QiyN) ■ (28) 

k=l 1=1 

Representing the 5-function in terms of a Fourier integral over a (see (jHl)), the integrations 
over yi . . .y^ can be easily performed: 

Wiy) = J ^ exp{i{y-x)} exp | - y '^afaf 



I 



Performing the a integration finally fields: 



Qiy) - ^ exp {- [ , a' = J2afaf . (29) 



I 

,2 



The optimal result is achieved if the error squared, i.e., a , is as small as possible. Here, 
we must take into account the normalisation condition in fl7r|) . We therefore minimise 

^afcrf — \ai — min , (30) 

I 

where A is a Lagrange multiplier. The global minimum is easily obtained: 

A 2 1 

"^ = 2^' A =4-^- ^^'^ 

The minimal value for cr^ then satisfies 

^" = ^ ^ A = EA- (32) 

2 ^ erf 

The optimal choice for the weights can therefore also be written as 

ai = ^ . (33) 

Let us return to the above example. We find 

a ^ 0.089 , (34) 

ai ^ 0.30 , aa ^ 0.12 , as ^ 0.19 , ag ^ 0.39 . 

With the weights at our disposal, we easily find the optimal value for the world average 
y ~ 3.059. Together with the error in ( 134|) . the final result is 

X = 3.059 ± 0.089 = 3.06(9) . (35) 

Note that the true result x = 3 is again covered within error bars and that the error became 
significantly smaller than that of the best result provided by the BNL group. 
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1.3 Autocorrelations 



In the previous subsections, we have repeatedly assumed that the measurements Xi are 
independent. We will see below that a vital tool of computational quantum field theory is 
to use information from the measurement Xi to obtain the value for Xj+i. In this case, the 
data set Xi, i — 1 . . .n is generated by the chain 



Xi > X2 ^ • • • ^ 3^n— 1 ^ -^m 

and the probability of finding a particular set does not factorise anymore: 

P{xi,...,Xn) ^ Pixi) ... P{Xn) . 

In the context of QFT simulations we will, however, make an effort to render the values 
Xi as independent as possible. This generically imphes that events which are separated by 
some 'timeV, i.e., the events Xi and be considered as statistically independent. 

The trick for obtaining an idea of the error of the estimator is to group b measurements 
together: 



1 



M 



i=l 



where we choose 



1 < r < 6 



(36) 



(37) 



This latter constraint implies that the values are statistically independent and that 
they have a Gaussian distribution because of the central limit theorem. The quantities of 
interest are the average 



y 



1 ^ 



(38) 



i/=i 



which converges to the observable x in the limit M — > (X), and the corresponding error 



1 



n 



(39) 



where the cumulant C2[Py\ is given by 



' ^ M 



i/=l 



(40) 



with 



/n 
Yl dXi f{xi ...Xn) P{xi, ...,Xn). 
1=1 



Assuming translational invariance, i.e.. 
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we find 

^ b b 

(^^[Py] = p XI XI ^ ^) ' C{k-i) = (XiXk) - {Xi){Xk) , (41) 

i=l k=l 

where c{k — i) is called the autocorrelation function. Introducing the relative distance 
t = k — i, and trading the summation over k in (HTl) for a summation over t, we find 



= (42) 



i=l t=l-i 



Interchanging the summation indices t and i and after summing over i, this last equation 
becomes: 



C2[Py\ = ^C(O) + |5^(6-t)c(t). 



i=l 

We have already mentioned that the measurements Xj and (almost) uncorrelated. 

The equivalent statement is that the correlation function vanishes for sufficiently large 
arguments: 

c{t) ^ , for t > r . 
For b ^ T, we approximately find: 

C2[Py] ^ ^C(O) + I '^bc{t) = i . (43) 

i=l t=l-b 

It is convenient to introduce the normalised autocorrelation function by 

Pit) = ^ , c(0) = (x^) - (x)^ = c,[P^] . (44) 



The integrated autocorrelation time is then defined by 

b-l , 6-1 

2 



r = ^ E /^('^) = ^ + • (45) 

t=l-b t=l 

Inserting fl45ll44|l43p into (!39l) . we finally obtain for the error which should be attributed to 
our estimate y in (l38l) : 

= ^ C2[PJ . (46) 

Let us perform a consistency check by considering the special case that the measurements 
are uncorrelated. In this case, the autocorrelation function p{t) vanishes for t > 1, and the 
autocorrelation time is given by r = 1/2. We indeed recover the familiar result 

0"^ = — C2[Px] , for independent measurements. (47) 
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Note that autocorrelations increase the error bars. Not knowing the autocorrelations in a 
numerical simulations leads us to the erroneous assumption that the error of the estimator is 
given by (HTll . while the true result is by a factor r larger. Not knowing the autocorrelations 
always leads to an underestimation of the statistical error. 



2 Lessons from the Ising model 
2.1 Phase transitions 

A phase transition occurs if the properties of matter change qualitatively when an external 
parameter, such as temperature, is altered. The phase transition of water from a liquid 
to a gas phase when the temperature exceeds roughly f» 100° Celsius (under normal 
conditions), is well known from everyday life. A second example is the ferromagnet: the 
interaction between microscopic spins favour a unique orientation of the spins. This yields 
an ordered phase at low temperatures. Above the critical temperature, called the Curie 
temperature in the present context, the spins are organised in a disordered phase. 

Let us assume that a ferromagnet is in the disordered phase at a temperature slightly 
bigger than the critical temperature T^. If we decrease the temperature, the information 
of the unique orientation spreads over the spin lattice. This 'Gedankenexperiment 'shows 
that the spatial correlation of the spins becomes large at the critical temperature. This 
phenomenon, is quantified with the help of the spin-spin correlation function: 



a{x) a{y) ) oc exp 



\x 



y\ 



(48) 



The correlation length ^ obviously measures the spatial distance over which the spins show 
roughly the same orientation. Close to the phase transition, i.e. for T > T^, ^ becomes 
large anticipating the ordered phase: 



(T > T, 



(49) 



where z/ is a positive number. The divergence of the correlation length at the phase 
transition is characteristic for a transition of 2nd (or higher) order. In the case of a 1st 
order transition, the increase of ^ is hindered by the nucleation of bubbles which contain 
chunks of the new state of matter. These bubbles provide additional disorder and the 
correlation length stays finite. 

For phase transitions above first order, the singularity of the correlation length has its 
fingerprint in many other thermodynamical quantities such as the specific heat C or the 
magnetic susceptibility x- 



C ^ Co 



1 - 



T 



X 



Xo 



T 



(T > T, 
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The critical exponents z/, a, 7 are independent of the microscopic properties of the spin 
model (such as the lattice geometry), and only depend on the symmetries (present at the 
transition) and the number of dimensions. They are often used to sort solid state physics 
models into the so-called universality classes. 



2.2 Quantum field theory rising 

Let the lattice spacing a denote the distance between two neighbouring lattice sites . In the 
previous subsection, we found that the correlation length diverges if the coupling constant 
/3 (or inverse temperature in the present context, (3 = 1/T) approaches its critical value 
(see ( H9|l ). This statement can be phrased in units of the lattice spacing as 

^ = «: (^/?c - , /3^/3c, (50) 

where k is a dimensionless constant which can be obtained by numerical means. 

Let us now reinterpret these findings. Rather than saying that C, diverges and a is fixed, we 
say that the correlation length ^ is fixed and is given by an observable in physical units. We 
will see that this interpretation of the same data defines a quantum field theory. For 
fixed correlation length ^, ( 150|) defines the lattice spacing as a function of P, i.e., a a{P), 

a{(3) = l(/5-/?,)-e, P^Pc (51) 

K 

The key point is if we make the number of spins bigger and bigger and, at the same time, 
the distance a between the spins smaller and smaller, we will obtain a field theory in the 
limit a —>■ 0. For the 2d Ising model on a cubic lattice, we have (3c ~ 0.44 and u = 1. The 
field theory limit is then approached when a vanishes linearly with (5 — f3c- 

a{p) = i (/5 - , (2d Ising model). 



Note that the dimensionless parameter (3 is not at our disposal anymore, since it specifies 
the magnitude of the lattice spacing. Instead, the value of ^ parameterises the emerging 
quantum field theory. The exchange of a dimensionless parameter for a scale dependent one 
in a quantum field theory is known as dimensional transmutation. It is a generic feature 
of quantum field theories. For instance in the case of perturbative QCD, the dimensionless 
gauge coupling g is eliminated in favour of the scale dependent parameter Aqcd- 

Let us assume that a certain correlation function was obtained by a numerical simulation 
of a classical lattice model for large values \x — y\, 

D(^\x-y\^ = ^F(0(x)) F(0(y))^ oc exp|-m |x - y|| , (52) 
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Figure 2: Spin correlation along the diagonal and the symmetry axis, respectively. 

where m is called the screening mass. Since the distance \x — y\ is only known in units 
of the lattice spacing by construction, the simulation will provide the mass in units of 
the lattice spacing as a function of /3, i.e. ma(/5). If universality holds, one finds the 
characteristic scaling of the lattice model, i.e., 

ma(/?) = - , /5 < . (53) 

Using (1501) . we see that the product approaches a constant in the vicinity of the critical 
limit: 

= ma— = K^nK . (54) 
a 

Note that k and are two c-numbers which we obtain from the numerical simulations. 
With the help of these two numbers we can "measure" the desired mass m in units of 1/^, 
where ^ is the only free parameter of our theory. 

In the case of a quantum field theory, we expect that due to the isotropy of the vacuum 
the correlation function ^% only depends on the distance between x and y. In the clas- 
sical lattice model, continuous rotational symmetry is violated due to the presence of the 
cubic lattice, and it might happen that the quantum field theory which emerges from the 
lattice model inherits the anisotropy. This anisotropy can be measured by comparing the 
correlation length in lattice units along a lattice symmetry axis, ^, and a long the diagonal 
direction, (see figure [2]). As far as global symmetries are concerned, the symmetry under 
consideration is restored in the critical limit (ISTl) : 

id = i , for a . 

Further details on the restoration of rotational symmetry in the context of the 2-dimensional 
Ising model can be found in [2]. 
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2.3 Mean- field approximation 

The starting point for a thermodynamical description of the Ising model is the partition 
function: 

Z = 5^exp(-/3i7(a)) (55) 

where (3 = 1/T and where a spin cr^ = ±1 is associated which each site x of the square 
lattice. The sum in fl3^ extends over all possible spin configurations. The ferromagnetic 
interaction favours a unique orientation of the spins, and is described by 

H{cr) = - ^ a^ay , (56) 

<xy> 

where the sum extends over all pairs < xy > of nearest neighbours. In order to preserve 
translational invariance, periodic boundary conditions are often used in particle physics 
applications, although these conditions are difficult to interpret in the solid state physics 
context. 

In order to gain an initial insight into the phase structure of the Ising model, we choose a 
particular spin of the lattice, and assume heuristically that that we might replace the 
neighbouring spins by the mean 

(a) = i 5^ a,, expi-pHia)) . (57) 

Wx} 

The Hamiltonian is then approximately given by 

H{axo) ~ const. - 4(cr) ■ (58) 

Note that each spin possesses 4 neighbours on a cubic 2d square lattice. Equation fl37|) 
turns into a self-consistency equation to determine the (a), which can be interpreted as 
the magnetisation per site: 

i^) = ^ J2 expi-pH{a)) , (59) 

X = Yl exp(-/3if(a)) . (60) 

Performing the sum over leaves us with a non-linear equation: 

(a) = tanh (4/5(a)) . (61) 

Before we proceed with a numerical solution of this equation, we point out that ( 16T1) always 
possesses the trivial solution 

(a) = . 
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p 

Figure 3: Magnetisation per site as function of the inverse temperature (3\ solid line: mean 
field approximation; dashed line: exact 



A graphical inspection of (]6T1) easily shows that for 



(62) 



two non-trivial solutions ±c, for c > exist. The physical interpretation of the solution 
is clear: for sufficiently small temperature (high j3), an ordered phase exists. The critical 
value is, in mean-field approximation, given by 



MF 



1/4. 



(63) 



Equation (lUTj) can be easily solved numerically with the Newton method or by fixed point 
iteration. The result for the magnetisation as a function of the inverse temperature is 
shown in figure [3l Also shown is the exact result [3,4]: 



1/8 



1 - 



sinh^(2/?). 



(64) 
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Figure 4: Geometrical objects on a lattice. 
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The mean-field result qualitatively reproduces the correct phase structure. The mean field 
approximation is able to describe the transition from the disordered to the ordered phase. 
However, the mean field approximation fails at a quantitative level. The correct value for 
the critical value, which was already obtained by Kramers and Wannier in 1941 [3], is given 
by 

f3c = 0.44068679... (65) 

is significantly underestimated. Also the rise of the magnetisation close to (3c is not correctly 
reproduced. A Taylor expansion of ( l6Ti) with respect to /? around = 1/4 (and therefore 
also with respect to a), yields: 



{cr) 



(66) 



The mean field critical exponent of 1/2 is much too large compared with the exact exponent 

of Kxact = 1/8. 

The advantage of the mean-field approximation is that it can be easily applied to a variety 
of models (e.g. the Ising model in ci > 2 where no exact results are available). It often 
provides a correct first impression of the phase structure. The disadvantage is that it is 
difficult to improve the approximation in a systematic way. 



2.4 Duality transformation 

Let us list different geometrical objects on a lattice. The sites on a lattice are labelled by 
integer coordinates. Links are short line segments which join two neighbouring sites on 
the lattice. In order to unambiguously address a link on the lattice, we use coordinates 
which are integers with the exception of one coordinate which is half integer, such as 2.5 
(see figure H] for an illustration). Another important object is the so-called plaquette , 
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which is an elementary square of the cubic lattice. Two coordinates are half integer when 
a plaquette is addressed. In higher dimensions, there are also cubes, and their coordinates 
are half integer, while the other coordinates are integer. 

The dual lattice is an important object which helps to gain non-perturbative information 
for certain lattice models. The coordinates of the dual lattice are obtained by adding 0.5 
to all coordinates of the lattice. If we consider a d dimensional lattice model, the duality 
transformation maps an x-dimensional geometrical object into a d — x dimensional object 
on the dual lattice. Let us consider 2 dimensions. A site, such as (2, 4) is mapped into 
(2.5, 4.5), which are the coordinates of a plaquette, while a link, e.g. (1.5, 5), maps into 
another link namely (2, 5.5). 

With these prerequisites, let us consider the probabihstic measure of the 2d Ising model. 
Since the product a^ay can only be ±1, we expand: 



a + b a^Cy . 



exp|/3cr^aj^| 

Inserting both possible values for the product (7x<Jy, we find: 

a + 6 = e^ , a - b = e'^ , 

and finally: 

exp|/9 (Ji^o-yl = cosh/? + sinh P (TxCTy . (67) 
Hence, the partition function in ( |55l) can be written as 



Z = ^cosh^^ n + tanh/5 a^ay . (6. 



where x and y are nearest neighbours on the lattice, and the corresponding link is denoted 
by {xy). Note also that, in 2 dimensions, there are 2N links for a lattice with N sites. 
In order to perform the sum over all spin configurations in ( !68l) . we use the important 
relations: 

(T=±l cr=±l 

Hence, if we perform the sum over the spin in flUSj) . we must make sure that it appears 
twice (or an even number of times) when we expand the products of the square brackets. 

Thus, if we avoid a vanishing contribution to the partition function, integration over the 
spins generates closed loops the corners of which are marked by a pair of spins. For each 
link of the closed loop, we get a factor tanh/5. Hence, after we have integrated out all 
spins, the partition function can be written as: 

N(L) 



Z = cosh^^l3 2^ ^[tanh/3 , (69) 



loops 
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Figure 5: Integration over spins generate closed loops. 



where N{L) is the number of links of the closed loop L. Note that we obtained a factor of 
2 for each sum over a particular spin a. This gives rise to the prefactor 2^ in front of the 
sum in ( |69ll . We now have converted the Ising model into a string theory, but we have not 
gained much information on the Ising model so far. To proceed further, we must control 
the sum over closed loops. For this purpose, we introduce new variables r = ±1 which are 
associated with the plaquettes (see figure E]) • If we consider two neighbouring plaquettes, 
there is always just one link between them. Now we say that if the product of the two 
neighbouring plaquettes is —1, the corresponding link is part of the loop. If the product 
is 1, the link is not part of the loop. The advantage of the r variables is that we can 
randomly assign ±1 to them and all loops which we produce are closed. Hence, summing 
over all possible r configurations will do the sum over all possible closed loops for us. 

Note that each plaquette of the lattice is mapped to a site on the dual lattice. The link 
between two neighbouring plaquettes is mapped into the link between the adjacent sites 
of the dual lattice. Finally, we must express N{L) in terms of the r variables. For this 
purpose, we have to count all activated links (links which are part of a loop) on the lattice. 
It is easy to check that 

{xdyd) 

counts these links: if two neighbouring rs are equal, they do not contribute to N{L), and 
is they are different, they contribute 1 as they should. Using the r-representation of the 
closed loops, the partition function fl69|) becomes 



Z = [cosh/?] 2^ [2 tanh/3]^ J] JJ ^tanh/3 



J '^x d '^y d 



(71) 
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Figure 6: Introducing dual variables to represent closed loops. 



This last equation can be written as 

Z = sinh^(2/3) exp{ ^ 

(72) 

p = - ^ In tanh p . (73) 

We have obtained again a 2d Ising model which is now formulated on the dual lattice: the 
only difference is that the coupling is now (3 rather than j3. It is not generally true that 
the duality transform yields the same lattice model just with different couplings. Models 
which do have this property are called self dual. 

Now let us assume that /3 is large (small temperature). In this case, we find from (I73|) that 

P ^ e"^^ , f3 large. 
By contrast, if /? is small (the high temperature limit), we find 

P ^ — ^ InP , P small. 

Hence, large P corresponds to small P and vice versa (see figure [7|). This is interesting since 
the so-called strong coupling expansion techniques are available for small p. Performing 
the expansion with respect to P in the dual model, the large P regime can also be studied 
by analytic methods. The basis of this expansion is a Taylor expansion of the exponential 
with respect to p. This expansion naturally reaches its radius of convergence when P 
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128x128, 2D Ising, Wolff cluster alg. 




Figure 7: The dual coupling constant /3 as a function of P (left). Magnetisation as a 
function of /? (right). 

approaches the critical coupling j3c- Performing the expansion using the dual model, the 
Taylor expansion with respect to (3 also breaks down at the critical coupling. There are no 
other couplings for which singularities in thermodynamical quantities occur. Hence, the 
critical point is obtained if 

/3 = /3 = /3c. (74) 

Using (1731) . we therefore find 

/5c = - ^ In tanh/3c , f^'= = \ ^^(1 + ^2) ~ 0.44068679 .... (75) 

Figure [7] also shows the magnetisation as a function of (3 for a 128 x 128 lattice compared 
with the mean field result and the exact result in the infinite volume limit. 

There are lot of interesting features of field theories already present in the Ising model: 
there is the relation between a lattice model and a theory of strings, and there is the duality 
transform which maps the high temperature theory onto a low temperature theory. 

3 Markov chain Monte-Carlo: the Ising case study 
3.1 Foundations 

The idea central to all simulations of lattice models is to generate lattice configurations 
{(Jx} according to their probabilistic measure 

Pia) = expi-pH{a))/Z (76) 
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where Z is the partition function (1551) . A straightforward idea to accomphsh this task would 
be to generate randomly the spins at each site x and to accept or reject this configuration 
according to (|76l) . The problem is that we would hardly find any acceptable configurations. 
Why is this so? 

Let us answer this question in the context of the Ising model of the previous section. The 
two dimensional lattice consists of = 125 x 125 sites. Since o G {— 1,+1}, there are 
2^ ^ ]^g4704 fiiffgj-gj^t lattice configurations. We further introduce the average action per 
site, i.e. 

1 / 



N 



where 



^ <xy> 

A{x) : = 



Aix] 



A 



(77) 
(78) 



y>x,\x—y\=l y>x,\x—y\=l 

and where we have used translational invariance. A measure for the strength of the flue 
tuations of the action around its average value A^ s is given by 



K<xy> 



Ns 



-\ 2\ 



5^((^(x) - A) {A{y) - A)) . 



(79) 
(80) 



x,y 



The crucial observation is that the connected correlation function 

D{x-y) := {{A{x) - A) {A{y) - A) 

exponentially decreases for large values of \x — y\, i.e. D{x) oc exp{— x/^yi}, where is 
the correlation length characteristic for fluctuations in the action density. Hence, one finds 
that its integrated strength, the so-called susceptibility p, is finite at least for j3 ^ /3c, 
i.e., 

p := E D{x) < oo . (82) 

X 

These findings tell us that the standard deviation 6 ( IHOj) linearly grows with the number 
of sites, i.e. 6'^ = N p. 

Using the central limit theorem to estimate the probability for accepting an action density 
s, we find 



exp 



{Ns - Nsf 
52 



exp 



is - -s) 

p 



N 



(83) 



Hence, in the case of many sites, only configurations with an action per site close to the 
average action density can significantly contribute to the partition function. If we randomly 
choose the spins on the sites the action density can take any value between —1 and 1, and 
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the argument (s — s)/p is generically of order 1. Hence the acceptance rate is down to 

g^l28xl28 10^'^^^^ 

The basic idea is to only generate configurations which are relevant. Starting from a seed 
configuration cq, we will generate subsequent configurations Ci, C2, . . . , where the result 
for Cj+i should only depend on the precessing configuration and must not depend on the 
configurations Q-i. In this case, the set of configurations, 

Co > Ci > C2 > C3 > . . . > Coo 

is called a Markov chain. Central ingredient to a Markov chain is the probability W{b, a) 
with which configuration b is created out of configuration a. This probability must satisfy 
certain constraints: 



(i) Normalisation J2b W{b,a) = 1, Va 

(ii) Ergodicity W{b,a) > 1, Va, 6 

(in) Stability J2a ^(^' ^(«) = ^i^), V6 , 

where P{a) is given in (1761) . If these conditions are met, the series q converges to a con- 
figuration which is distributed according to P(coo) (EEj). In order to see this, we introduce 
the probability Qi{c) for finding a configuration c at position i of the Markov chain, and 
denote the deviation from the desired distribution by 

e^ = 5Z|Qi(c) - P(c)| . (84) 

c 

Because of property (ii), there is a Wmm with 

W{a,b) > W^in > , W'{a,b) := W{a,b) - W^in > . (85) 
Furthermore, the condition (i) implies that 

5^Q,(c) = 1, as 5^P(c) = 1. (86) 

c c 

Using the stability condition (iii), we then obtain: 

c a c a 

= J2\j2^V'{c,a) b,(o) - P{a)] + - P{a)] 



c, a 



Qi{a) 



P(a) 



^7) 
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Using the triangle inequality and the positivity of W, we find 

em < 5^5^1^'(c,a) Q,{a) - P{a) 



Changing the order of summation and using (see ( l85l) ) 

^iy'(c,a) = J2^W{c,a)-W^,^) = 1 - n,o„f 1^^;, , (89) 

c c 

where riconf is the number of configurations, we finally find convergence: 

Ci+l < [1 - '^confW^min] ^ |(5i(a) - P{a) = [1 - Wconf W^min] ■ (90) 

a 

Instead of demanding the less stringent condition (iii), one often demands detailed balance: 

(iii)' W{b,a) P{a) = W{a,b) P{b) . 

The latter condition immediately leads to condition (iii) if we sum the equation (iii) ' over 
the configurations a: 

Y,W{b,a)P{a) = J]W^(a,6)P(6) = P{b) , 

a a 

where we have used condition (i). Since condition (iii) follows from (iii)' and only (iii) is 
necessary for our proof above, demanding detailed balance, i.e., (iii)', is more restrictive. 

3.2 Heat-bath algorithm 

The heat-bath algorithm works as follows: (i) randomly choose a site Xq and consider the 
corresponding spin o"(xo) for the update. Since the spin only interacts with its nearest 
neighbours, the interaction can be written as 

H = const. - hoCTxo , ho = ^ a^c ■ (91) 

<xxo> 

The relative probability for choosing a{xo) = 1 is given by exp{hol3}, and the relative 
probability for cr(so) = —1 is given by exp{— /iq/^}. (ii) Calculate the absolute probability 

exp{-/3ho} + exp{/3h(,} 

with which the spin must be set to 1. Choose a random number z G [0, 1]. For z < p, 
choose = 1 otherwise choose = —1. (iii) Subsequently, pick another spin for the 
update and start again with (i). Once all spins have been visited at least once, one sweep 
has been performed. 
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Figure 8: Thermalised spin configurations of tlie 2d Ising model for several (3 values 
starting from high temperature phase to the low temperature ordered phase. 



The algorithm above needs an initial configuration. We could choose a unique orienta- 
tion of all spins. Since this is the ground state of the Hamiltonian which dominates the 
partition function for small values of the temperature, this initialisation is called a cold 
start. Alternatively, we could start with a random orientation of the spins. This is a con- 
figuration which is relevant at very high temperatures where interactions are negligible. 
This initialisation is therefore called a hot start. Independently of our choice, a number 
of sweeps is carried out to generate a statistically relevant configuration. This procedure 
is known as thermalisation. The number required to arrive at an equilibrated spin lattice 
depends on the number of degrees of freedom and the temperature. 

Let us now examine typical lattice spin configurations. Starting at low the sample 
configurations are highly disordered (see figure [H]). Increasing /5 up to ~ 0.3, the clusters of 
spins with the same orientation already extend over several lattice spacings. Approaching 
the critical value, e.g. for (5 ~ 0.42, the clusters are already as large as the lattice. This 



25 



observation reflects the growth of the spin correlation length which, for the present case, 
is 

-1 

(T > T,) , (93) 
and hence diverges when /3 — /?c- 

If our numerical approach should produce the configurations of a Markov chain, the con- 
figurations may not depend on the Monte-Carlo history. To find out whether the config- 
urations are indeed statistically independent, we may inspect the autocorrelation time r 
e.g. say for the magnetisation M (recall subsection II. 3l for discussions of autocorrelations). 
To guarantee independence, we perform about 2r Monte-Carlo sweeps before we consider 
a configuration eligible for contributing to an estimator. In the case of the heat bath algo- 
rithm (in fact, for all local update algorithms), one discovers that the autocorrelation time 
strongly increases when the critical point is approached. This implies that the interesting 
regime of the model, namely the regime close to the phase transition, is not accessible 
with these types of algorithm. The reason for this is the following: consider a spin inside 
one of the clusters. All the neighbouring spins are pointing in the same direction. If this 
spin is now subjected to a local update procedure, the spin hardly changes because of the 
strong mean field produced by the other spins. Hence, only the boundaries of the clus- 
ter are significantly modified after one sweep through the lattice. The correct physics is, 
however, described by configurations consisting of strongly fluctuating clusters. In order 
to change a cluster completely, there are roughly lattice sweeps necessary. Hence, only 
after sweeps, the configuration has changed significantly. This, however, implies that 
the autocorrelation time is roughly given by r ^ Indeed, it was empirically observed 
for the Metropolis algorithm that 

r ^ r , ^Metro- 2.125. (94) 

The index z is called the dynamical critical exponent and depends on the algorithm. Since 
the physical correlation length ^ diverges at the phase transition, flMl) implies that the 
regime near the phase transition cannot be simulated with local update algorithms. 



T 



3.3 Cluster update algorithms 

State-of-the-art simulations which explore the transition regime use the so-called cluster 
algorithms. The difference to local update algorithms is that many spins are flipped at a 
time. To derive the prescription of such a cluster update, we rewrite the partition function 
dSSD as 

Z = ^exp I /3 ^ a^ay | = XI 11 CFxCTy) . (95) 

{(Tx} \ <xy> / {cr^} <xy> 

If both neighbouring spins, and ay, in ( 195|) are equal, the probabilistic factor in ( 195|) 
equals exp f]. For an opposite orientation of the spins, the probabilistic factor is given by 
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exp(— /?). We therefore cast into 



^=En 



1 -p) + p5<,,a. 



p := 1 — e' 



-2/3 



(96) 



We now are going to make the representation of the partition function even more involved 
by using the identity 



a + 6 = 



n=0 



a5no + bSni 



Introducing variables rixy € {0, 1}, which are associated with each link of the lattice, we 
obtain 



^=EE n 



=/3 



Wx} {nxy} <xy> 



(97) 



The cluster update prescription is now obtained by performing standard heat bath steps 
for the variables {(Tx) and {^xj/}- 

Let us consider the update for the link variables Uxy first. In order to avoid generating a 
configuration of vanishing probability, we must choose n^y = if the neighbouring spins 
{(Tx) and {ay} are different (see (jSZ]))- If these spins are equally oriented, the probabilistic 
measure in (P7|) ) is given by 

implying that the link Uxy is set to 1 with probability p. Given an initial spin distribution, 
the values of all link variables can be chosen according to the above prescription. 

Let us now consider the spin update. According to probabilistic measure, i.e., 

~l~ V ^(JxO-y ^nxy,l J 

all spins which are connected by links nxy = 1 must be of equal orientation. All spins 
which are connected by so-called activated links, said to be part of a 

cluster. The task is now to find all such clusters on the lattice. Once these clusters have 
been identified, we assign ±1 (with equal probability) to all spins of the same cluster. 

This first versions of such cluster update algorithms are due to Fortuin and Kasteleyn [5], 
Swendsen and Wang [6] and Wolff [7]. It is found empirically that the dynamical critical 
exponent is strongly reduced: 



T ~ 1^ , ^cluster ~ 0.2 

Introductory discussions can be found in [5-8]. 



(9J 
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4 Quantum field theories on computers 



4.1 Quantum mechanics 

Let us assume that the motion of a particle of mass m in 1 dimension is governed by a 
potential V{x). The classical equation of motion can be calculated by variational methods 
from the action ^ 

S = j^dt l^x^ - V{x)^ . (99) 

Classically, these equations of motion determine the time evolution of the position of the 
particle x{t). At quantum mechanical level, the partition function 

Z{T) = tr exp|-^i7| (100) 

is a convenient starting point to discuss the thermodynamics of the physical system. Here, 
H is the quantum mechanical Hamiltonian, i.e.. 



T is the temperature, and is considered as an external parameter. Once one has succeeded 
in calculating the partition function fllOOp . thermodynamical quantities can be easily ob- 
tained by taking derivatives, e.g., the temperature dependence of the internal energy is 
given by 

, , ^9 d In Z(T) 

{H) = T' . (102) 

Although a direct calculation of the eigenstates (n| of the Hamiltonian might be the easiest 
way to calculate a quantum mechanical partition function in practical applications, I would 
like to reformulate fllOOp in terms of a functional integral. This will be of great help when 
we generalise the quantum mechanical considerations to the case of a quantum field theory. 

For these purposes, I introduce a length scale L := 1/T and an interval [0, L] which I 
decompose into equidistant portions of length a <^ L, where a is called lattice spacing. 
It is trivial to obtain 

expj-^i/l = expl-^ a h\ = JJexp{-ai/}. (103) 

Let us define complete sets of momentum and position eigenstates {\p) and respectively) 
by 

1 = dX„ \X„) {X,y\ , I = I dpy \Py) {p ^\ , (l04) 



ioi u = 1 . . .N . As usual, these states obey 

{Pk\xk) = expj-^pfcXfcj 
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Using a complete set |xo) of position eigenstates to evaluate the trace in fllOOl) . we find 



/ dxo {xo\ JJ^ exp{— ai/jixo) = / dxodpodxidpi ... dxjy-idpN-i 



{xo\e " \pq) {po\xi) {xi\e " \pi) {pi\x2) ■■■ (xAr_i|e " \pn-i) {pn-i\xq) . 
Note that the operators p^ and V{x) do not commute. We may, however, write: 



exp<j-a- aV{x) + —\y{x),p^\ + 



p 

exp {— al^fa;)} exp <( —a 

2m 



Since \x) and \p) are eigenstates of the position operator and momentum operators, re- 
spectively, we find 



(xfcl exp{-aif}|pfe) = exp < -a 



2 

1^ + V{xk) + 0{a) 
2m 



exp{-pkXk} 
n 



The partition function therefore becomes up to terms of order 
Z(T) = / dxo (ipo dxi dpi . . . dx^-i dp^-i dxjyexp < —a 

I fc=0 

exp <^ -^^Pk{xk+i - Xk) > (a;o|a;Ar) 



2m 



(105) 



k=0 



It is straightforward to perform the momentum integrations, which are Gaussian, 
Z(T) = ( — — j / dxo dxi... dxN ^ox^v 



(106) 



N-l r 



exp <^ -a 



k=0 



m { xk+i -Xkf 
a? h? 



+ V{xk) 



This equation is a completely regularised expression for the partition function and can 
be directly used in numerical simulations. Note that in the framework of quantum field 
theory, one sets h= 1. 

A compact notation can be derived by formally taking the lattice spacing a to zero. For 
this purpose, we define ah := ha, and the Euclidean action by 



m 



dT{ '-^X^ + V{X) 



(107) 



Note the sign change in front of the potential compared with the standard action ( l99i) . 
The interval [0,L], which was introduced above fll03p . is called Euclidean time interval. 
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Figure 9: Classical versus quantum partition functions of a 1- dimensional particle chain. 



By construction (see above), the length of the Euclidean time interval is given by the 
inverse temperature, i.e., L = 1/T. We also introduce a Euclidean particle trajectory, and 
a Euclidean velocity 

x,~^x{t) ^m^^i(r), (108) 

ah 

where we identify dr = ah- Using the shorthand notation 

I / dxo dxi . . . dxM^i — > Vx(t) , 

\ CLh J J 

the partition function fll06p can be formally written as a functional integral 

(109) 



Z(T) = j Vx{t) exp|-^5£;| . 



Eq. fll09p suggests that an average over all Euclidean trajectories x(r) must be performed 
where the probabilistic weight of each trajectory is given by exp{— S'^;/^}. Note also that 
in view of the 5-function in (11061) only trajectories which are periodic in Euclidean time 
must be considered, i.e., a;(0) = x{L = 1/T). 



4.2 Quantum field theory 

For illustration purposes, we consider the 1-dimensional particle chain in figure M Here, 
the positions of the particles i = . . .n are characterised by their extensions Ui from the 
equilibrium position. The particles experience a harmonic potential depending on the 
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distance to the nearest neighbour. Here, I choose the boundary conditions Mq = 0, 
The Hamiltonian, which describes the classical physics, is given by 



0. 



n-l 

E 

1=1 



2m 



Pi 



;iio) 



111) 



Hence the classical partition function is given by the multi-dimensional integral 

In order to calculate the full quantum mechanical partition function of the particle chain, 
we first write down the Euclidean partition function. Note for this purpose that the 
displacements now acquire an additional dependence on the Euclidean time Ui — Uiij) = 
Uri- With this notation the Euclidean action is given by 



N n-l 

EE' 

T=l 1 = 1 



m 
2^ 



+ 



D 



[u 



'112) 



The interactions between the c-number fields Uri can be easily visualised (see figure [9]): the 
fields Uri harmonically interact with their nearest neighbours. The harmonic interaction 
strength is given by D/2 in space direction and m/2a^ for neighbours in Euclidean time 
direction. The quantum mechanical partition function can be calculated by integrating 
over of the fields Uri located at the sites of a 2-dimensional grid, .i.e. 



Z{T) oc j Vu exp{— S^;} 



;ii3) 



where the temperature enters the consideration via the extension of the lattice in Euclidean 
time direction with fields obeying periodic boundary conditions. 

To conclude, we observe that the partition function of a classical D + 1 dimensional field 
theory (in lattice regularisation) describes the full partition function of a D dimensional 
quantum system. D is the number of space dimensions. This correspondence is very 
helpful in understanding the quantum behaviour of a theory, since it can be mapped to a 
classical field theory (at the expense of an additional dimension). In the next section, we 
will study classical partition functions in 4-dimensional Euclidean space in order to derive 
the information on the thermodynamics of the full quantum system. 



5 Lattice gauge theory 
5.1 The gauged Ising model 

The Ising model, strictly speaking the partition function fl55p . is invariant under the global 
transformation of the spins given by 

a^{x) = n a{x) , n = ±1. (114) 
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The transformation is called global because the transformation affects all spins at the same 
time, i.e., Q is independent of the coordinates (sites). The corresponding symmetry group 
is Z2. 

This symmetry group can be generalised to a local symmetry, also known as gauge sym- 
metry, by demanding invariance under 

a^{x) = n{x) a{x) , n{x) = ±1 . (115) 

Of course, the action (1561) of the standard Ising model is not invariant under the huge 
symmetry group which is now [^'2]^, where is the number of sites. In order to obtain 
a version of the Ising model which possesses a Z2 gauge symmetry, we need to change the 
action. The only way to do it, is to introduce an additional field, Z^{x). This field is 
associated with the links of the lattice: x specifies the site and n the direction in which we 
find the link. Alternatively, we could write: 

Zf,{x) = Zi^^y) , y = X + , 

where is the unit vector in [i direction. For the latter expression, we will also abbreviate 

a; + = x \ \i . 

For the action, we choose 

•S'matter = X] ^{^) Zp{x) o{x + /i) , (116) 

and demand that the link Z^ transforms under gauge transformations as 

z^{x) = n{x) z^{x)n{x + fi) . (117) 

Since spin and link transform simultaneously with the same Q{x) and since Q'^{x) = 1, one 
easily proves the gauge invariance of the action flll6p . 

Obviously, the action ^matter describes the interaction between the matter fields, i.e., the 
spins, and the new link fields. What is left to do is to design a gauge invariant action for 
these new degrees of freedom. This interaction should be short ranged in order to preserve 
some desirable features such as universality. A possible choice is 

5'link = Pfiuix) , P^,u{x) = Z^{x) Zy{x + fl) Z^{x + v) Zy{x) . (118) 

Here, the numbers (x, n > u) specify the plaquette of the lattice the lower left corner 
of which is located at site x and which is spanned by the directions and v. The field 
combination P^u{x) is often called the plaquette for short. The proof that Pfj.u{x) is indeed 
invariant under gauge transformations flll7p is left to the reader. 
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The total action of the gauged Ising model consists of two parts: the matter part and the 
"gauge" part. Correspondingly, there are two coupling constants: the convention is that (3 
is the pre-factor in the pure gauge action, while k multiplies the matter part. 

Once our system is now gauged, it only makes sense to consider gauge invariant observables 
since non-gauge invariant quantities vanish. Let us explore this for a simple gauge variant 
quantity such as the spin correlation function: 

C{Xo,yo) = ^ (^{Xo)(7{yo) exp ^S[a]^ , S[cr] = ^matter + ^'unk , (119) 

Z = ^exp{s[a]}, (120) 
{^} 

Let us now consider a particular gauge transformation flllSp of the spins, i.e.. 

Renaming all spins in the sum in (I119p by cr(x) — cr^(x), we use the gauge invariance of 
the action and the sum , i.e., 

S[a] = S[a^] , E = E • 

{a} {a^} 

The sum is trivially invariant, since we sum anyhow over all possible ±1 combinations for 
the spins. Thus, we obtain: 

C{xo,yo) = ^ E ^"(a;o)cr"(l/o) explsfcr]! 

= \ E e^pj'^M} = -C{x^.y^), (122) 

where we used our particular choice 01211) in the last line above. We conclude from this 
that C{xq, yo) = 0. 



5.2 Pure Z2 gauge theory: 3 and 4 dimensions 



Let us now consider the particular case k = when the matter fields are absent from the 
theory. The emerging theory is called pure Z2 gauge theory, and it is a theory of link fields 
only. What are the physical (i.e. gauge invariant) degrees of freedom in this case? Let us 
consider the more interesting case of 3 and 4 dimensions for these considerations. In order 
to talk about gauge invariant information, we now consider the plaquettes Pp, p = (x; /iz/), 
defined in (11181) . We say that a short flux line (vortex) passes through the plaquette p if 
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Figure 10: Flux passing through an elementary cube (left) and closed flux lines on a 3d 
lattice (right). 



Pp = —1. Since the plaquette variables P are gauge invariant, so are the flux lines. More 
formally, we introduce a vortex plaquette variable by 

Vp = l[Zi (123) 

and consider the flux lines which enter/leave an elementary cube of the lattice. We take 
the product of all vortex plaquettes which are associated with the faces of the elementary 
cube and find 

pGc 

where u is the total number of vortices at the faces of the cube. Inserting the definition 
( 11231) . we also find that 

pGc pGc iGp 

since in the latter products all Zi factors appear twice (see figure [TOl left panel; remember 
Zf = 1). Comparing f ll25p with (11241) . we realise that u must be even. In particular, u = 1 
is excluded implying that a vortex never ends inside a cube. Considering 3 dimensions 
(or the spatial hypercube of 4-dimensional space time), we find that the gauge invariant 
vortices form closed lines in space. See figure [TUl right panel, for an illustration. 

4 dimensions: Let us consider the 4 dimensional model first. The constraint (11251) is 
most easily interpreted on the dual lattice. A plaquette p of the original lattice maps onto 
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a plaquette *p of the dual lattice, and a cube c corresponds to a link */ on the dual lattice. 
Hence, the constraint (11251) reads on the dual lattice 



l[v,p = 1 (126) 

*pe*z 

This simply means that the number of vortex plaquettes which are attached to a link on 
the dual lattice must be even. Accordingly, the vortex plaquettes form closed surfaces on 
the dual lattice. If n denotes the number of negative plaquettes on the original lattice, the 
number of trivial plaquettes in 4 dimension is 6A^ — n, where is the number of lattice 
points. Hence, the probabilistic weight of such a configuration is: 

[exp {Z?}]'"^-" [exp{-/5}]" = exp{6N(3} exp{-2/3n} , 

so that the partition function can be written as 

Z = exp{6N(3} exp{-2(3n} . (127) 

{closed surfaces} 

Let us interpret this partition function: the degrees of freedom are closed two dimensional 
sheets (2-branes) embedded in four dimensions. The surface A of these branes is given by 
n. Hence, the probabilistic factor is given by 



exp 



-2(3 a} 



implying that 2f3 can be interpreted as the surface tension. At zero temperature {(3 oo), 
the empty vacuum (no 2-branes) is realised. At finite temperatures, the brane entropy 
competes with the penalty from the weight factor. A direct calculation of the entropy 
of 2d world-sheets in 4d would be cumbersome. However, exploiting the relation to the 
Z2 gauge theory makes the calculation of brane expectation values easily accessible by 
numerical means. 

Let us proceed to obtain the duality map of the 4d Z2 gauge theory. The basic trick 
to perform the sum over the closed surfaces is to introduce degrees of freedom which 
automatically resolve the constraint. In the present case, these are links *Z^:i on the dual 
lattice. Let us consider 

5^ + tP*p[*Z]] , (128) 

where P^p[*Z] is the plaquette generated by the links *Z on the dual lattice. When we 
remove the brackets in fll28p . the only way to have a non- vanishing contribution to the 
sum is by making sure that each link *Z^:i appears an even number of times. This, however, 
means that the negative plaquettes P^p[*Z] form closed surfaces. Hence, we find 

5Z n[i + ^p*prz]] = 2'"^ (129) 

*P {closed surfaces} 
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Note that there are 4N hnks on the 4 dimensional lattice and that the factor 2^^ arises 
from the sum over *Z^i. Thus, using fll29p in (11271) . we find 



Z = exp{6N(3} 2-^^ H + tP*p[*Z]\ , t = exp{-2/?} . 

Since the plaquette P^p[*Z] only acquires values ±1, we may write: 

JJexp|/3P,p[*Z]} = [cosh^]^^ n + tanh^P,p[*Z] 
*p *p 

The partition function (ll3Up therefore becomes 

2 ^ r exp{/j} l^^ 
. cosh P . 

exp{— 2/3} = tanh/3. 



(130) 



(131) 



2-^^ E exp{/3 5^P.,rZ]}, 



(132) 
(133) 



First of all we note that the dual of pure Z2 gauge theory is another 4-dimensional Z2 
gauge theory: the model is self-dual. Furthermore, relation fll33p is already familiar to us: 
we have obtained a complete analogue of the relation between f3 and its dual f3 for the 2d 
Ising model. We therefore once again encounter the fact that the weak coupling regime is 
mapped to the strong coupling regime of the same model. As a byproduct we find that 
the critical coupling is given by 



= I ln(l + V2) 



0.44068679 . . . 



(134) 



The fact that the critical couplings of the 2d Ising model and 4d pure Z2 gauge theory 
coincide might be a numerical accident. At least, I do not know any deeper reason for 
it. I finally point out that the Z2 gauge symmetries of the original and the dual formu- 
lation are completely unrelated. This can be most easily seen from the fact that, at an 
intermediate stage, we have formulated the model entirely in terms of physical, i.e., gauge 
invariant variables: the closed vortex sheets of the dual lattice. Resolving this constraint, 
the gauge invariance of the dual model arose from a parameterisation invariance, namely, 
the redundancy when we performed the sum over all closed world sheets with the help of 
dual gauge links *Z. 

3 dimensions: Let us finally discuss the 3-dimensional model. In 3 dimensions, a 
plaquette p is mapped to a link */ and a cube c is mapped to a site *x on the dual lattice. 
The constraint f ll25p then translates to 



n 



v*i 



(135) 
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The plaquettes carrying negative flux on the original lattice are represented by links forming 
closed loops on the dual lattice. The partition function now takes the form (there are 3N 
links on the lattice): 



Z = exp{3NP} J2 exp{-2/3n} 
{closed loops} 



(136) 



We already know how to perform the sum over all closed loops from subsection EI 
2-3N Y\ exp|/3r,,r,j,} (137) 

_n 3N 



Yi n 

cosh pi 3JV 



^ JJ^l + tanh/5 

{t,x} *l 



cosh P 



J2 [tanh/5]" 
{closed loops} 



Identifying once again 
we find 



exp{-2/3} = tanh/? 



Z 



exp{/3} 
.2 cosh^. 



3N 



expj^ Pn^ny^ . 



(138) 



(139) 



Obviously, the Z2 gauge theory is dual to a theory which is not a gauge theory anymore, the 
3d Ising model. This has tremendous consequences: for the standard Ising model, cluster 
update algorithms are available. Using the Swendsen-Wang or Wolff type cluster update, 
we are able to simulate a gauge theory with much less autocorrelations. Unfortunately, 
such a framework is not (yet) available for more relevant theories such as lattice Yang-Mills 
theories. 



5.3 Setting up lattice Yang-Mills theory 

Due to the universality conjecture, the lattice model with the correct number of dimensions 
and the correct symmetries uniquely defines the corresponding quantum field theory in the 
critical limit. The purpose of the present subsection is to propose a classical lattice model 
which satisfies these prerequisites in the case of Yang-Mills theory. 

The QCD matter fields (quarks) belong to the fundamental representation of the so-called 
SU{Nc) colour group {Nc = 3 for QCD). Gauge invar iance means that the action of the 
quark fields is invariant under the local unitary transformations, i.e., 

q{x) q'{x) = n{x)q{x) , n{x) G SU{N,). . (140) 

As explained in many text books, local gauge invariance of the quark kinetic term may 
only be achieved by introducing additional dynamical fields, the gluon fields Afj,{x). 
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Figure 11: Path ordered product of link variables. 



The quark fields are associated with the sites in a lattice formulation. Hence, the symmetry 
group of the classical lattice Yang-Mills model is [SU{Nc)]'^'' , where Ng is the number of 
lattice sites. In order to enforce such a high symmetry in the critical limit of a lattice 
model, it has turned out essential to realise the symmetry even for finite values of the 
lattice spacing a. This in turn forces the model to maintain local gauge invariance in the 
continuum limit [9] . A potential candidate for a quark kinetic term in the non-interacting 



case IS 



E 



(141) 



where the 7^ are the Euclidean 7 matrices. Of course, the action (11411) is not invariant under 
the gauge transformations (I140p . To achieve this invariance, we introduce an additional 
field of vector type, thus being related to the links of the lattice. 



(142) 
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plaquette 

Figure 12: Lattice plaquette variable 



Generalising the quark kinetic term (11411) to 

= ^\ Qix)lf^U^ix)q{x + fi) -g(x + /i)7^[/^(x)g(x) 



x,ii 



(143) 



one obtains the desired local invariance upon demanding that the link fields transform as 

Uf,{x) n{x)U^{x)n'^{x + 12) . (144) 

Let us follow the case of the gauged Ising model and construct a kinetic term for the link 
fields U^{x). For lattice models "kinetic" means that the interactions of the fields on the 
lattice are short range, i.e., only nearest neighbours are involved. In order to design a 
gauge invariant kinetic term for every value of the lattice spacing, we firstly investigate the 
transformation properties of a path ordered product of link variables. Let us consider an 
open path C which starts at point x and ends at y (see figure [TT] for an illustration), and 
define 

P{x,y) = ^ n ^(^) ' (145) 

where V implies path ordering. Inserting the gauge transformed links (11441) into (I145p . one 
finds 

P{x,y) ^ P'{x,y) = n{x)P{x,y)n{y) . (146) 

With the help of (I145p . it is easy to construct a kinetic term for the link variables which 
(i) is gauge invariant and (ii) involves only next to nearest neighbours. For this purpose 
one chooses C to be a closed path starting at x and ending aX y = x which encircles an 
elementary plaquette (see figure [T2|) : 



1 

1 

1^. 



tr P{x, y) 

tr \u^{x) U,{x + fi) UUx + u) Ul{x) 



(147) 
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Using fll46p and the invariance of the trace under cychc permutations, one easily shows 
that the plaquette (11471) is indeed gauge invariant. 

The lattice partition function involves an integration over the dynamical fields of the theory. 
In the case of the group valued link variables fll42p . the question arises which measure VUfj, 
should be employed for the integrations. We must demand that the integration measure 
does not spoil gauge invariance. To ensure this we use the so-called Haar measure which 
satisfies 

dUf,{x) = d(^AU^{x)B^ , A,Be SU{N^) . (148) 

The Haar measure is available in closed form for the unitary groups SU{Nc). Here, I will 
only present the Haar measure for SU(2) group integrations. The SU(2) unitary matrix U 
is conveniently parameterised in terms of the Pauli matrices, 

U = ao + iar, UU^ = 1 + = 1. (149) 

Since the constraint = 1, i.e. + = 1, is not changed if U is multiplied with A 
from the left and B from the right, respectively, these multiplications can be viewed as 
rotations in the 4-dimensional space spanned by (ao,a). Therefore, an invariant measure 
can be defined by 

dU = dao dai da2 da-s (51 Oq + — 1 I . (150) 



Introducing polar coordinates for the 3-dimensional vector a := an, nn = 1, the integration 
over the norm of the vector a can be performed with the help of the 6 function in (11501) . 
We obtain the final result for the SU(2) Haar measure, i.e.. 



dU = dao y 1 — dQ^ , (151) 
which is commonly used in lattice simulations. 

Finally, the lattice representation of the gauge invariant partition function is given by 

Z = jvUVq Vq^ expj-^g + [^/^-(^) + ^'^J | ' (1^2) 

where the quark interaction is encoded in Sq (11431) and P^u{x) is the plaquette (I147p . /3 
is related to the bare gauge coupling constant g of the continuum formulation hj P = 
2Nc/g'^. The particular choice (I152p of the lattice regularised gluonic action is known as 
the Wilson action [9]. Note that the fields q{x), g^(x) are anti-commuting Grassmann 
fields. This choice for the fermionic fields is necessary to obtain the correct Fermi statistics 
as well as to ensure the Pauli principle. It implies that the lattice model (11520 cannot be 
straightforwardly be used in numerical simulations. Rather, since the action for the quark 
fields is quadratic, the integration over the quark fields has to be performed analytically: 

j VqVq^ expj-g^M^sgB} = BetM[U] . (153) 
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where the index A comprises space-time as well as spinorial, etc. indices. The quark deter- 
minant DetM[f/] is a gauge invariant function of the link variables U^{x). Note, however, 
that the link interaction mediated by the quark determinant is non-local, implying that 
a link at a particular site is coupled to all other links of the lattice. In practice, this 
implies that a local update of a single link enforces the calculation of a functional deter- 
minant. This explains why the numerical simulation of Yang-Mills theory with dynamical 
quarks requires much more computational resources than the simulation of the theory in 
quenched approximation, where the quark determinant is neglected for the update of the 
link variables. 



5.4 The fermion doubling problem 

It turns out that the treatment of the quark degrees of freedom in (11521) is still too naive: 
since the Dirac equation is linear in momentum, its lattice analogue does not produce just 
the desired quark degrees of freedom in the limit a — 0, but rather 2^ copies of them {D 
is the number of space time dimensions). This is already true for the free theory as will be 
shown in what follows. 

Let us firstly introduce the generating functional for connected Green functions in the case 
of free and massless bosonic theory, 

%l^/l'^exp{-i,.n.,,^,.,.}. ,154) 

A sum is understood over indices which appear twice. One easily verifies that the connected 
correlation function is obtained from Z[j] via 

/(.-.):=(^.^.)-(^.)(^,)=^. (155) 



By "completing the square" in (11541) . we find 

Z\]] oc exp|^j^(^n"^ 



(156) 



and hence for the free bosonic case 

'n-M . (157) 



In order to evaluate the inverse 11 ^, of the "kinetic" operator, we introduce its eigenvalues 
and eigenvectors, whereupon 

Ii\k) = \k\k) , (158) 

and formally write 

(n-i) = E 1^) f (^1 • (159) 
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k a ka 



Figure 13: Dispersion relation from the tree level kinetic term (continuum versus lattice 
formulation) for the bosonic case (left) and the fermionic case (right panel). 



It is now easy to calculate the correlation function for the continuum case 11 = —9^. The 
eigenfunctions are subjected to periodic boundary conditions (j){x) = 0(a; + L), i.e., 



[x] oc e 



Akx 



AkL 



27r 

k = —n , n E Z . 



(160) 



The discrete k levels are called Matsubara frequencies. In the continuum, there are no 
further restrictions on the integer n. Making the ansatz fll60p . we find that the eigenvalues 
are given by 

\{k) = k'^ (continuum) . (161) 

Hence, a free massless particle manifests itself in the correlation function f ll59p as a pole 
at zero momentum transfer. The lattice version of the eigenvalue equation is 



-0(x + /i) + 2 0(x) - 0(x-/i) 



(162) 



In order to solve this equation, we use the plane wave ansatz f ll60p . One crucial difference 
between the lattice and the continuum version is that only wavelengths / obeying 



2 ^ " 



71 



(163) 



are sensible. The lattice naturally provides an UV momentum cutoff, i.e., Auv = ^r/a. 
Inserting (11600 into (11620 one finds 



(164) 



For momenta which are small compared to the UV cutoff, i.e., ka ^ vr, we recover the 
continuum dispersion relation 



A 



latt 



2„2\ 



1 + 0{k'a 



(165) 
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In figure [T^ the continuum dispersion relation for bosons is compared to its lattice version. 
The lattice correlation function has only one singularity reflecting that in the scaling limit 
\a? <^ 1, fca ^ TT, the dispersion relation of the continuum free particle is recovered. 

Let us move on to the fermionic case. In order to reproduce the correct Fermi statistics, 
fermion fields ip{x) are of Grassmann type and obey anti-periodic boundary conditions. I 
refer to the textbook [2] for an introduction to the free fermionic theory, and only quote the 
final result for the correlation function which formally agrees with f ll59p . In the continuum, 
the eigenvalue equation is given by 



(166) 



where anti-hermitian (Euclidean) 7 matrices are used. The ansatz for the spinor wave 
functions is again of plane wave type. 



ipi^x) oc u{k) e 



ikx 



AkL 



27r / 1 
-1 , k = —[n + - 



neZ 



The spectrum \{k) is determined by making the ansatz 



u{k) 



+ A 



Mo 



which yields 



and therefore 



iljt — X u{k) = ilj: — X ilj: + X 



k'-X' 



Uo 







(167) 



(168) 



(169) 



(170) 



Hence, the spectrum of the continuum theory is linearly increasing, A = itvfc^. Using 
the kinetic energy for a free quark theory introduced in fll4ip . the lattice version of the 
eigenvalue equation is given by 



X a ip{x) 



(171) 



The ansatz fll67p also provides the eigenvectors of the eigenvalue problem fll7ip . Repeating 
the steps which have led to the continuum dispersion relation, one finds its lattice analogue 



Xa = sin^ f fc^a 



(172) 



The fermionic eigenvalue distribution is shown in figure ffT^ . right panel. Close to the 
critical limit {Xa -C 1), one recovers the continuum dispersion relation from f ll72p by 
making a Taylor expansion with respect to ka. In addition, a second singularity occurs for 
ka ^ IT. This shows that even if Aa <^ 1 a second fermion flavour arises from the lattice 
fermion action (11410 . 
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It can be shown that this fermion doubhng problem must occur for a chirally invariant 
action which is translationally invariant and local (Nielsen-Ninomiya No-Go theorem). 
At the present stage, a lot of research effort is devoted to incorporate chiral symmetry at 
the expense of, say, a moderate non-locality of the action [10]. 



5.5 Overlap fermions 

In the continuum formulation, the chirally invariant Dirac operator V satisfies the relations 

{©,75} = 0, {V-\^^} = 0, (173) 

which tells us that the non-zero eigenvalues A appear in pairs {A, —A}. Let D denote a 
lattice candidate of the Dirac operator (in units of the lattice spacing a) satisfying the 
so-called Ginsparg- Wilson relation [11], 

{^,75} = 2D75D. (174) 

One observes that the right hand side of fll74p is of order (compared with the order a of 
the left hand side) implying that the naive continuum limit a — > of (11741) reduces to the 
chiral relation (11731) . The most important observation, however, is that a certain remnant 
of the chiral symmetry is still present in the lattice version. Defining 

D'^ := - 1 , (175) 

and using 

{75,D-^} = 275, 

which directly follows from the Ginsparg- Wilson relation (11741) . we observe that may 
be used as a chirally invariant quark propagator, i.e.. 

Hence, we are left with the task to find an operator D obeying the Ginsparg- Wilson relation 
(11741) . Here, I will briefly discuss the Overlap Dirac operator [12]- [14], firstly introduced 
in the pioneering paper [12] by Narayanan and Neuberger. One introduces 



D = i 

2 



1 + 



;i77) 



where if is a Hermitian operator with eigenvalues ±1. Common choice is 

1/2 



H = D^l {dID^^ , (178) 



where is the standard Hermitian Wilson-Dirac operator. Inserting (I177p into (11741) . it 
is straightforward to prove that D from (11771) satisfies the Ginsparg- Wilson relation (11740 . 
A comprehensive discussion of the quark propagator (I175p in the context of a simulation 
of SU(3) Yang-Mills theory can be found in [15]. 
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Figure 14: Wilson loop and loop-loop correlation function 



5.6 Measuring observables 

We have observed that the trace of the path ordered product P{x, y) of link variables f ll45p 
is gauge invariant when taken along a closed curve C, i.e., x = y. Depending on the choice 
for the loop C, the expectation value of such loop variables can be connected to physical 
observables. For instance, for the so-called Wilson loop, we choose a rectangular loop with 
size r in one spatial direction and the extension t in the Euclidean time direction (see figure 
[HI left panel). In the limit of large t, the Wilson loop expectation value is related to the 
potential V{r) between a static quark and a static anti-quark which are separated by the 
distance r, i.e., 

(w[C]'j oc exp|-V(r) t| , (179) 

In the particular case that the potential is linearly rising, V{r) = ar with string tension 
cr, one observes that the Wilson loop expectation value exponentially decreases with the 
area A enclosed by the loop C. Since a linearly rising quark anti-quark potential implies 
confinement (see discussion below), this area law (due to Wilson) is a litmus test for quark 
confinement. 

Furthermore, one can calculate the correlation function L{tx—ty, x — y) of two loops centred 
at X and y, respectively (see figure \\M . right panel). Here, information is transported from 
point X to ?/ by gauge invariant states \ph) . The shape of a particular loop determines its 
behaviour under the symmetry transformations of the underlying lattice. These symmetry 
transformations correspond to rotations in the continuum limit. Therefore, it is possible 
to select the spin quantum number of the state \ph) by adjusting the shape of the loop. 
For large distances A = t^, — ty, the correlation function exponentially decreases, i.e., 

L(ta; — ty, u = f — y) oc exp<— maAi. (180) 
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Figure 15: The static quark anti-quark potential as obtained from pure SU(2) lattice 
gauge theory. Plot from [16]. 

Hence, the calculation of loop correlation functions provides access to the so-called screening 
masses m of physical particles. In the purely gluonic theory, the only gauge invariant states 
are the glue balls, while in full QCD also hadronic states contribute to the correlation func- 
tions. 



5.7 The continuum limit 

For definiteness, I confine myself to the case of pure (i.e. no quarks) SU(2) gauge theory. 
The generalisation of the findings of the present section to SU {Nc) is straightforward. The 
task is now to find the critical limit of the lattice Yang-Mills theory. 

There is a lesson to learn from continuum Yang-Mills theory. In order to renormalise the 
continuum theory, one absorbs a logarithmic divergence into the bare gauge coupling. A 
detailed calculation yields 



11 
24^ 



A2 

In — + finite 

/i2 



:i8ii 



where A is the UV cutoff and where fi is an arbitrary renormalisation point. The coefficient 
in front of the logarithmic term is the quantity of interest and can be obtained by evaluating 
a bunch of one- loop Feynman diagrams. Eq. (11811) shows that in the critical limit A 
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Figure 16: Approaching the continuum hmit of SU(2) (left) and SU(3) (right) lattice 
gauge theory (improved action from [16]). 



oo the bare coupling vanishes. This is one manifestation of the celebrated property of 
asymptotic freedom. Switching from the continuum to the lattice formulation we identify 
A = vr/a and use P = 4/ g"^ to straightforwardly derive 

a^(P) = const. expJ-^/^l . (182) 

Due to asymptotic freedom, we expect that the critical limit is approached when /5 — > oo. 
The perturbative relation between a and (3 in (11821) is called asymptotic scaling. 

Modern computer simulations use a more complicated "kinetic" term for the gluon fields. 
One example of such an improved action is given by 



where PjiJ{x) is the 2x2 Wilson loop. Imposing the constraint 

+ 16 /t2 = 1 , (184) 

ensures that the familiar relation between (3 and the bare gauge coupling g, P = 2Nc/g'^, 
is maintained. The residual freedom of choosing ki and k,2 can be used to obtain a rather 
good agreement with asymptotic scaling on rather coarse lattices. 

In order to search for the critical limit with the help of numerical simulations, we calculate 
a physical quantity, e.g. the string tension a in units of the lattice spacing as a function 
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of the only parameter (3. This is done by calculating the static quark anti-quark potential 
V{r) as a function of the quark anti-quark distance r = na. The outcome in units of the 
lattice spacing is shown in figure [151 By fitting the numerical data to 



V{r)a = Vq — — + aa? n , 



n 



we find the string tension in units of the lattice spacing, aa^, for each value of (3. The 
outcome of this calculation is shown in figure [THl One indeed observes that the c-number 
aa? exponentially decreases for large values of (3 in agreement with the prediction (11821) of 
continuum Yang-Mills theory. The quantum field theoretical limit of the classical lattice 
model is obtained by interpreting the correlation length, i.e., the string tension a in the 
present example, as a fixed physical quantity, and reinterpreting the (3 dependence of the 
numerical data for ao? as the (3 dependence of the lattice spacing. 

Let us assume we have obtained a glue ball mass m in lattice units, i.e., we know ma as 
a function of (3. If the mass m is a physical observable, one must recover from the data 
the characteristic dependence a{[3) (see (I182p ) for sufficiently large [3 values. Hence, the 
ratio of the two dimensionless numbers m?a? /aa^ approaches a constant for (3 close to the 
critical point (see figure [TU| right panel). Extrapolating the data to the continuum limit 
a — > 0, i.e., (3 —>■ oo, one determines the physical mass m in units of another physical 
scale, i.e., ^/o^. Finally, let us count the number of parameters. The only parameter of 
the classical lattice model is (3, but (3 is no longer at our disposal in the quantum field 
theory limit (which implies f3 — >• oo). However, the physical value of the correlation length 
(or s/a in the present example) takes over the role of a free parameter. The replacement 
of a dimensionless parameter by a mass scale in the continuum limit is a feature of many 
quantum field theories and is called dimensional transmutation. On the lattice every 
mass scale is obtained in units of the string tension, y/a = 440 MeV is used to assign the 
familiar units of QCD to observables. For 32 lattice points in any space-time direction, we 
then find: 



For a fixed number of lattice points, we note that we cannot make f3 arbitrarily small since 
the physical volume becomes too small. Small values of f3 result in large volumes, but we 
cannot make (3 too small in order to have a reasonably large UV cutoff. Thus, for a fixed 
number of points, there is a small window of f3 values which are appropriate for a study of 
QCD particle properties. This window is sometimes called the scaling window. 

Acknowledgements: I thank Tom Heinzl and Martin Lavelle for a careful reading of 
the manuscript and helpful comments. 



P (input) 



1.250 1.400 1.500 1.600 



aa^ (calculated) 
L = Na 
A = vr/a 



0.279(2) 0.0922(7) 0.0528(3) 0.0311(2) 
7.7 fm 4.4 fm 3.3 fm 2.6 fm 
2.6 GeV 4.6 GeV 6.0 GeV 7.8 GeV 
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